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Many nanostructures today are low-dimensional and flimsy, and therefore get easily distorted. 
Distortion-induced symmetry-breaking makes conventional, translation-periodic simulations invalid, 
which has triggered developments for new methods. Revised periodic boundary conditions (RPBC) 
is a simple method that enables simulations of complex material distortions, either classically or 
quantum-mechanically. The mathematical details of this easy-to-implement approach, however, have 
not been discussed before. Therefore, in this paper we summarize the underlying theory, present the 
practical details of RPBC, especially related to a non-orthogonal tight-binding formulation, discuss 
selected features, electrostatics in particular, and suggest some examples of usage. We hope this 
article to give more insight to RPBC, to help and inspire new software implementations capable of 
exploring the physics and chemistry of distorted nanomaterials. 
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I. MATERIAL SYMMETRIES BEYOND 
TRANSLATIONS 

Translational symmetry and Bloch's theorem has been 
the backbone of materials research for a long time 1 . 
Bloch's theorem was originally associated with simu- 
lations of bulk crystals and translational symmetry in 
three dimensions, and this association is still strong. 
Together with periodic boundary conditions (PBC), or 
Born- von Karman boundary conditions, Bloch's theorem 
has enabled simulating the infinite bulk using a single, 
minimal unit cell. 

Nanoscience, however, has introduced novel material 
structures that often lack the translational symmetry. 
These structures include nanotubes, nanowires, ribbons, 
nanopeapods, DNA, polymers, proteins, to mention only 
a few examples. Some structures, such as nanotubes, 
are often modeled with translational symmetry, which 
is, however, dangerous because of their low-dimensional 
character and flimsiness — in experiments the real struc- 
tures get distorted and Bloch's theorem and conventional 
PBC becomes invalid. To remedy this problem, during 
the course of time different research groups have inde- 
pendently developed new methodological improvements. 

The seminal ideas to use generalized symmetries were 
presented by White, Robertson and Mintmire, who used 
helical (roto-translational) symmetry to simulate model 
carbon nanotubes (CNTs) using a two-atom unit cell 2 , 
which enabled calculating CNT properties that were 
inaccessible before 3 ' 4 . A similar approach, also ap- 
plied to CNTs, was followed by Popov 5 ' 6 . Liu and 
Ding used rotational symmetry to simulate the elec- 
tronic structure of carbon nanotori within a tight-binding 
model, although with a static geometry 7 . Dumitrica 
and James presented their clever objective molecular 
dynamics (OMD), using classical formulation 8 ; later, 
OMD was extended to tight-binding language 9 . Also Cai 
et al. presented boundary conditions for twisting and 
bending, using nanowires as their target application 10 . 
These approaches towards generalized symmetries have 
been applied to various materials and systems, such as 
Si nanowires 9 , MoS2-nanotubes 11 , nanoribbons 12 , the 
bending of single-walled 7 ' 13 ' 14 and multi-walled 15 CNTs, 
and vibrational properties of single- walled CNTs 16 . 



Recently, in Ref. 17, we introduced revised periodic 
boundary conditions (RPBC), a unified method to simu- 
late materials with versatile distortions. This method has 
a particularly simple formulation, with both classical and 
fully quantum- mechanical treatments. It can be used in 
conjunction with molecular dynamics and Monte Carlo 
simulation schemes, it is designed for general distortions 
and all material systems, and it can be regarded as a 
generalization of previous work 2 10 . Because we only pre- 
sented the overall outline of RPBCs in Ref. 17, the pur- 
pose of this paper is, above all, to be the mathematical 
and technical companion of that work. We give in-depth 
formulation of the approach, suggest few examples of us- 
age, and discuss the long-range electrostatic interactions 
carefully, as to provide a detailed-enough overview to aid 
implementing and applying RPBC in practice. Exam- 
ples of earlier simulations with RPBC include bending 
of single- walled CNTs 18 , twisting of graphene nanorib- 
bons 17 ' 19 , and spherical wrapping of graphene 20 . 



The idea of using general symmetries in material simu- 
lations should not be surprising, as in physics and chem- 
istry symmetry has always played a special role. What is 
surprising, though, is that symmetries beyond the trans- 
lational ones have not quite reached the mainstream 
of materials modeling. Using generalized symmetries 
does not imply material distortions; they can be used 
also merely to reduce computational costs. For exam- 
ple, RPBC is based on revised Bloch's theorem, which 
in turn is based on the validity of Bloch's theorem for 
any cyclic group — an age-old group-theoretical common 
knowledge 21 . However, it is not until nanoscience and its 
low-dimensional distorted structures that have brought 
the motivation to finally go beyond the conventional 
translational symmetry, especially in practical simula- 
tions. We hope this paper could demonstrate how simple 
the RPBC formulation is, and what kind of ingredients 
are required in the implementation. Using symmetries 
beyond translation should be in reach also for the simu- 
lation community mainstream. 
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II. REVISED BLOCH'S THEOREM 

Let us first, to have a common starting point, repeat 
the translational Bloch's theorem, as it is found in all 
solid state textbooks 22 . In its conventional form, Bloch's 
theorem 

^k(r-T n )=e^ k - T ^ak(r) (1) 

is valid for a system of electrons in a periodic potential 
V(r) = V(r + T n ), where T n = mTi + n 2 T 2 + n 3 T 3 
and Ti are lattice vectors. In Eq. (1), the wave function 
^ak(r), labeled with band index a and k- vector, is the 
solution to Schrodinger equation, 

ffyak(r) = [p 2 /2m + V(r)] ^ ak (r) = ^ak(r), (2) 

where V(r) could also be the Kohn-Sham potential of 
the density functional theory. The system is subjected 
to periodic boundary conditions (or Born-von Karman 
boundary conditions), which mean that translation by 
Jjj = MjTj, the system's length in direction j, leaves 
the system unchanged. What has made the theorem (1) 
so powerful, is that the knowledge of the wave function 
within a single unit cell is sufficient to solve the electronic 
structure and structural properties of the crystal as a 
whole. 

Next we restate Bloch's theorem (1) in a form easier to 
generalize. First, we introduce the notation T n r = r + 
T n for the coordinate transformation of translation, with 
the inverse transformation T n r = r + T n = r — T n . 
Hence the left-hand side of Eq. (1) becomes V ; ak(7" _n r). 

Second, we define an operator D(T n ) acting on wave 
functions that induces the real-space translation T n . For 
this, recall that shifting a function "forward" [action by 

l)(T n )] is equivalent to shifting the coordinate system 
"backward" (V = r~ n r), 

Z)(r n )^(r)=^(T- n r). (3) 

The way to arrive at (3) is the following. Let us make 
a forward translation T n from the state i/j in the coor- 
dinate system r to a state ijj' in the coordinate system 
y' = T n r. The wave function does not change, meaning 
— ip'ir') or i/j(r) = ?//(T n r), which holds for any r. 
In particular, it holds for r = T _n r // , so (dropping the 
double prime) ifj(T~ n r) = ?//(r). Now, since the state 
is the state induced by the forward translation, we define 
D(T n ) by D(T n )^(r) = ^'(r), which leads to (3). With 
this notation the Bloch's theorem becomes 

£(T>ak(r) = MT- n r) = e" ik T " Vak(r). (4) 

Now we proceed to the revised Bloch's theorem, as it 
was presented in Ref. 17. It holds for a system of elec- 
trons in any external potential invariant under general 
(isometric) symmetry operation r' = <S n r, that is, when 

V(S n r) = V(r) (5) 

for any set of commuting symmetry operations <S z n % 
<S n = S^S^ 2 • • • . Using <S n instead of T n in Eq. (4), the 



revised version of Bloch's theorem, valid for any symme- 
try, reads 

D(S n )i> aK (r) = ip aK (S- n r) = e~ iKn Vw(r). (6) 

Here, the vector k generalizes the conventional k- vector 
and the vector of integers n = (m, n 2 , . . . ) gives the num- 
ber of transformations <S n = S^S^ 2 • • • for each symme- 
try Si. Knowing wave function Vw( r ) f° r arbitrarily 
chosen simulation cell is sufficient to calculate the wave 
function for atoms belonging to any (n th ) image, and 
thereby to solve, again, the electronic structure of the 
whole system. 

Let us next clarify the conditions when the revised 
Bloch's theorem (6) holds. Because the system should 
be left unchanged, the symmetry operations should com- 
mute; therefore, the operators D(Sj J ) must commute as 
well, 

[D(S;n,D(Sl*)} = 0. (7) 

Condition (5) is equivalent to [V(r), D(S n )] = and 

since D(S n ) commutes also with p 2 , requirement of po- 
tential invariance (5) becomes 

[H,D(S n )}=0. (8) 

The periodic boundary conditions are generalized in 
the following way. First, write the condition ^ak(r — 
Lj) = ^ak( r ) in terms of the translation operator, 
D{T M ')^{v) = ^{T~ M 'v) = Vak(r), where Mj = 
Lj/Tj is the number of lattice points along direction j. 
(Lj is the the length of the entire crystal in direction 
j.) Next, replace the symmetry operation T — » S and 
k y k. This gives £)(<S M )Vw(r) = Vw(r), or 

D(S M ) = D(sf [ S™ 2 . . . ) = 1, (9) 

where in general M = (M-[, M2, . . . ), with Mj being 
or Mj, where Mj is the number of transformations upon 
which the system is mapped onto itself. (HjMj = J\f is 
the total number of unit cells in the crystal.) 

In mathematical terms, Bloch's theorem deals with 
symmetries of the Hamiltonian alone and follows from 
the group representation theory for cyclic groups. In- 
deed, considering, for simplicity, the one-dimensional 
case, n —y n, the set of transformations 

{1, DiS 1 ), D(S 2 ), . . . , D(S M ~ 1 ) I D(S M ) = 1} (10) 

forms a cyclic group with one-dimensional representa- 
tion e i27rm / M ,ra = 0,l,...,M-l. Due to condition (8) , 
this group is also a symmetry group of the Hamiltonian. 
Therefore, the eigenfunctions of the Hamiltonian trans- 
form according to representation of its symmetry group, 
which is equivalent to Eq. (6). 

Since any unitary operator can be uniquely repre- 
sented by an exponent of a hermitian operator 23 , we 

have D(S) = e~ ik . Then, D(S n ) = e~ ikn , or in multi- 
dimensional notation 

D(S n ) = e~ ik n . (11) 



The operator k is commonly called the generator of the 
group. Due to condition (8), k commutes with H, and its 
components kj commute among themselves. Therefore, 
eigenfunctions of the operator with k\n) = k,\k), form 
a common set of eigenstates with the Hamiltonian opera- 
tor. The vector k, is the eigenvalue of the operator k and 
thus a good quantum number that can be used to label 
the energy eigenstates [as evident already in Eq.(6)]. The 
physical meaning of vector k relies on the symmetry of 
the system and is thus specific to given symmetry (see 
Subsection II A). 

Now, consider the eigenvalue problem for £)(<S n ), 

D(S n M)=e- ik - n \iP)=C n \*P), (12) 

where C n is a constant. We close Eq. (12) with (k\ to 
find 

(k\ e~ ikn = e~ iKn ^(k) = C„^(/s), (13) 

yielding C n = e~ lKn with i/j(k,) = (k^). Next we note 

that it is equivalent to define D(S n ) by ^(r) = ^'(r') = 

D(S n )ilj(S n r) and by D(S n )\r) = |<S n r); the unitar- 

ity property i)+(<S n ) = D'^S 11 ) = D(S~ n ) results in 

(r\D(S n ) = (<S -n r|. Taking this into account, we close 
Eq. (12) with (r| to find 

(r\D(S*m = <<S- n i#> = e— n (r|^), (14) 
or in coordinate representation 

D(S n )^j(r) = ^j(S- n r) = e~ iK n ^(r). (15) 

Applying the cyclic condition (9), we arrive at set of al- 
lowed values for k\ e tKjM ^ = 1 = e l27rrnj or 

/ mi m 2 \ . . 

n = 2tt [ — — , — -, . . . , (16) 

where m 3 =0,1,..., Mj — 1. 

Finally, labeling the eigenfunction with band index 
a and the good quantum number Eq. (15) becomes 
Eq. (6). 

We thus proved the revised Bloch's theorem of Eq. (6): 
for a system of electrons in the external potential result- 
ing from any symmetric arrangement of atoms (or any 
other external potential for that matter), wave functions 
separated by <S n transformations differ only by a phase 
factor. Therefore, by calculating the wave functions in 
a single unit cell — whatever its form — , one can simu- 
late the electronic structure of the symmetric system as 
a whole. 

Alternatively, the revised Bloch's theorem can be for- 
mulated by saying that the Hamiltonian eigenfunctions 
are of the form 

^a K (r)=e-- n « UaK (r), (17) 

where u aK (S m r) = u aK/ (r) is a periodic function and n(r) 
is a generalized dimensionless coordinate, chosen in a way 
to satisfy 

n(<S m r) = n(r) + m, (18) 
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FIG. 1: (Color online) Revised Bloch's theorem: (a) The 
product of the modulating phase factor, e* k r , and a periodic 
function u a \^(r) gives the translational Bloch state, (b) The 
very same principle applies to other symmetries. 



with m = (mi, m2, . . . ) being a vector of integers, so that 
(17) with condition (18) would satisfy (6) by construc- 
tion. The choice of n(r), however, is not unique. For 
most transformation n(r) can be a continuous function, 
yielding a fraction of the unit transformation for given r; 
for example, in one-dimensional case (n — >• n), at the ori- 
gin, which also gives the beginning of the simulation cell, 
n(r) = 0, at the middle of the simulation cell n(r) = 0.5, 
at border between the simulation cell and the first copy 
of the simulation cell n(r) = 1. With improper transfor- 
mation works as well, but requires a discontinuous n(r). 

For translations in Cartesian coordinates, we have 
rij(r) = Xj/Tj and Kj = 27rrrij/Mj [Eq. (16)], which 
yields k • n(r) = ^ 27rrrij / (MjTj) Xj = k • r, and we 
recover the well-known Bloch state 

^ak(r)=e ik - r Wak (r), (19) 

where i£ a k(r) = ^ak( r + T n ) is a lattice-periodic func- 
tion. General idea behind formulation (17) is illustrated 
in Figure 1. Next, we illustrate this formulation using 
selected familiar symmetries. 



A. Illustrations with familiar symmetries 

A simple transformation to consider is a two-fold ro- 
tation (7r-radian rotation). Since twice repeated it yields 
the original system (that is, S 2 = 1, M = 2, and 
m = 0,1), the ft-point sampling is n — and n = tt. 
Therefore, 

D(S n )i> aK (r) = 4> aK (S- n r) = (±l)"^ OK (r); (20) 
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the wave function for k = it becomes negative upon odd 
number of transformations. 

For M-fold rotation around z-axis in cylindrical co- 
ordinates r = (p,ip,z), S rn (p 1 (p 1 z) = (p,ip + ma,z), 
where a = 2tt/M is the angle of rotation. The condi- 
tion (18) is satisfied with n(r) = cp/a, and we further 
take k — 2i\mjM = am to write Eq. (17) in terms of az- 
imuth a 1 angle ip and quantum number m (which is noth- 
ing but the familiar magnetic quantum number): 

VWn(p, ^ Z ) = e im<P U am (p, Lp, Z), (21) 

where u am (p, (p + a n ,z) = u arn (p, cp, z) for any a n = na 
(n = 0, ±1, . . .). For example, with a = 7r, we arrive at 
Eq. (20). By considering a six- fold rotation symmetry, 
we can simulate the 12-atom benzene CqHq by 2-atom 
generalized unit cell (one CH unit); we know the wave 
function of the whole system if we know the wave func- 
tion ip arn (p,(p, z) in the wedge-shaped region p G [0, a], 
a = 7r/3. On the other hand, the spatial region for simu- 
lation cell does not need to be connected; in benzene, for 
example, the CH unit could be constructed from spatial 
region for C and H even on opposite sides of the molecule. 

A helical transformation with pitch length L z = MT Z1 
is a combination of a rotation by an angle x = 2tt/M 
and a translation by T z where the translation goes along 
the axis of rotation (here 2- axis), S m (p,p, z) = (p, p + 
ma,z + mT z ). The condition (18) can be satisfied with 
different choices for n(r), for example, n(r) = z/T z , 
n(r) = ip/Xi or n ( r ) = ( Z /T z + y?/x)/2, meaning that 
z/T Zl LpjXi or their combination can serve as the dimen- 
sionless coordinate. Further, take k = 2-KvnjM = mx to 
write Eq. (17) as 

Vw(p, <p, z) = e ™(W(Mr l)+ ,)/2 Uam{pj ^ ^ (22) 

where u arn obeys symmetry via the relation u arn (p,(p + 
a n , z + T n ) = u am (p, <p, z) for any T n = nT z and a n = nx 
(n = 0,±l,...). 

The above examples of proper symmetry transforma- 
tions can be described in terms of the corresponding 
generators k (Table I). The foundations of the revised 
Bloch's theorem, as this section has shown, are very fa- 
miliar. 



III. THE APPOINT SAMPLING 

One of the central practical issues in the RPBC ap- 
proach is the sampling of the k- values (^-points). It is 
essentially similar to the familiar k-sampling with trans- 
lational symmetry, although there are some differences. 
Depending on symmetry, the values of a given compo- 
nent Kj of the ^-vector may or may not be quantized, 
and sampling of kj hence falls in two schemes, in either 
discrete or continuous sampling. 

First, in discrete sampling the component 
Kj accepts only the values given by Eq. (16), 
{2imij /Mj \rrij = 0, 1, . . . , Mj — 1}, where Kj corre- 
spond to symmetry transformation Sj. This usually 
means a small Mj, say Mj = 2 or Mj = 6, but even if 
Mj would be large, say thousand, all thousand values 
do not necessarily need to be sampled. Non- allowed 
ftj's result in nonphysical wave functions and ultimately 



problems in numerical evaluations. Discrete sampling 
suggests that the corresponding symmetry is genuinely 
periodic, and the periodic boundary condition is a real 
physical condition. The simplest example of this case is 
the two- fold transformation discussed around Eq. (20). 

Second, in continuous sampling Mj goes to infinity (or 
can be treated that way), such that k,j = 27rmj /Ap- 
points can be sampled freely between [0, 27r), or within 
the Brillouin zone [— 7r,7r). Continuous sampling hap- 
pens for symmetry operations involving translation, for 
which the periodic boundary is not a real physical condi- 
tion, but rather a convenient mathematical trick that has 
turned out useful. Note that in regular three-dimensional 
system PBC means all dimensions to be periodic in an in- 
tertwined and bogus fashion; in two-dimensional system 
PBC represents topologically a toroid. 

Since RPBC works also with translational symmetry, 
there must be a relation between the ^-vector and the 
conventional k- vector; this was already shown in Ref. 17. 
It is obtained by equating the exponential factors in 
Eqs. (4) and (6), e" ik Tn = e _i ^' n . For given k, by 
choosing n such that only rij = 1, we get 
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k • n = Kj = Tj • k = T jA U = !> 2 > 3 ); ( 23 ) 
1=1 

here I runs though Cartesian components of vectors Tj 
and k. Solving Eq.(23) hence yields ^(k) or vice versa. 
It is easy to check that plugging k = Y^=i m j /Mjhj 
into Eq. (23) yields the condition (16), where bj are the 
reciprocal lattice vectors with bj • Tj/ = 2n5jj>. Eq. (23) 
suggests also that the ^-points are sampled continuously 
when k-point are. 

Usually for two symmetries, say Si and Sj, the sam- 
pling schemes are independent, meaning that, for exam- 
ple, K{ can have discrete sampling while Kj has continuous 
sampling. Yet sometimes the symmetry transformations 
can be coupled, which leads to coupling of K{ and Kj sam- 
plings. This can happen when system is such that map- 
ping of simulation cell atoms l\ times with «Si is identical 
to mapping of simulation cell atoms I2 times with 52, or 
S 1 ^ = S l 2 , which leads to the coupling of ^-points for 
symmetries S\ and 52- From the theorem (6) it follows 
that e _mi/l = e _m2 ^ 2 , and the coupling of the sampling 
is given by the condition 

l\K\ = I2K2 + 2i\m (m integer). (24) 

Here values of K\ are governed by Eq. (24) where m = 
0, 1, . . . , h — 1 and K2 acts as independent parameter 
[alternatively, K\ may act as an independent parame- 
ter; then, values of K2 would be given by Eq. (24) with 
m — 0, 1, . . . , I2 — 1]. The sampling of the independent 
parameter K2 [or K\] follows either of the two sampling 
schemes discussed above. In the general case we have 
S ni = 5" J , and the couplings of the sampling schemes 
can be obtained from e~ zn ' ni = e -^ n i+ 27rm . 

In the next section, we list cookbook-type recipes how 
certain symmetries can be used with the RPBC ap- 
proach. For each symmetry we point out the pertinent 
K-point sampling scheme. 
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TABLE I: Generators of one-parameter symmetry transformations. Translations are generated by momentum operator pj, 
rotations by angular momentum operator L z , and helical transformations along z-axis by p z and L z together. 
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IV. SELECTED SYMMETRY SETUPS 

The symmetry of the system determines the coordinate 
transformations <S n and the shape of the simulation cell. 
Every particle in the simulation cell, located at R/, has 
M images located at R" = <S n R/ with n = (ni, n 2 , . . .). 
A fairly general symmetry transformation for particles' 
coordinates is the affine transformation 

l S n r = ^(w n )r + T n , (25) 

where T n = niTi + n 2 T 2 + ^3T3 is translation and 
7Z(u> n ) is rotation for an angle uj n = niXi + ^2X2 + ^3<^ 
around an axis given by the vector u> n . This is not, how- 
ever, the most general form [cf. Eqs. (34)-(38)]. When 
rotation is a part of transformation, there can be at most 
one translation that should be collinear with the rotation 
axis (u? n || Ti). At the same time, the number of transla- 
tional transformations along this direction (possibly for 
different lengths) is unlimited; so is the number of coaxial 
rotational and helical transformations. Selecting some of 
the angles xi? X2, & and vectors Ti, T2, T3 yields the 
special cases we discuss next (see also Table II). The 
different symmetries are referred to, from practical sim- 
ulation viewpoint, as simulation setups. The titles in the 
following refer to the nomenclature of these setups. 



A. Bravais: translational symmetry revisited 

Transformations of this setup consists of up to three 
translations [xi 5 X2 5 ^ = in Eq. (25)]: 

<S^r = r + n,T„ i = 1,2,3; (26) 

S n = S^S^Sz 3 . This is the typical triclinic setup, 
mostly used for solids with Bravais lattices. The Ap- 
points, as the k-points, can be freely sampled. [The di- 
rect relation between a given K-point and a k-point is 
given by Eq. (23)]. 

Note that, for demonstrative purposes only, one- 
dimensional translation can also be constructed from two 
symmetry operations of the form (26), with Ti = Li and 
T2 = 2Li, <S n = S^S^ 2 - In this case the ^-points are 
coupled due to the relation S% = £2- Eq. (24) holds, 
yielding K\ = ft 2 /2 + 7rm, Mi = 2, and M 2 — >• 00. This 
kind of treatment is, of course, silly and artificial, but 
demonstrates the flexibility of the RPBC approach. 



B. Wedge: cylindrical symmetry 

Main transformation of this setup is rotation 
[T2,T3,xi,X2 = in Eq. (25); notations of Si are kept 
consistent with the discussion following Eq. (25)], 

S^ 3 r = n(n 3 cx)r, (27) 

possibly combined with one translation, 

<S x ni r = r + mTi; (28) 

<S n = S^S^ 3 . Translational periodicity in axial direction 
is either present (a || Ti 7^ 0, Mi —too), or not present 
(Ti = 0). This setup is for cylindrically-symmetric sys- 
tems. The parameter a is the angle of the wedge-shaped 
simulation cell. For a large angle a, k^'s are strictly dis- 
crete (a = 27r/small integer), while for sufficiently small 
angle a, k^'s can be sampled freely (the smallness is 
somewhat subjective, but a < 27r/20 is "small" for most 
practical calculations). 



C. Chiral: helical symmetry 

Transformation of this setup is helical [a, X2> T 2 , T 3 = 
OinEq. (25)], 

S^E^^r + niT, (29) 

where rotation and translation are coupled (xi II Ti). 
This setup is for chiral systems. The parameter xi is 
an angle of shift between two closest chiral units at dis- 
tance T z apart. The helix makes a full turn in 2tt/xi 
transformations and its pitch length is (2ir/xi)T z . k,i 
can be sampled freely, since the system is infinitely long 
(Mi 00). 



D. Double chiral: helical symmetry 

Helical setup can be extended with second helical 
transformation [a, T 3 = in Eq. (25)]: 

STr = ft(niXi)r + niTi, (30) 
S 2 n2 r = ft(n 2 x 2 )r + n 2 T 2 ; (31) 

both angles xi? X2 and both Ti, T 2 can be different (xi II 
% 2 || Ti || T 2 ). This setup is for chiral systems that 
do not need the whole circumferential part to represent 
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TABLE II: (Color online) Application of RPBC for distorting materials or improving computational efficiency. Some familiar 
and conventional examples, not particularly special for RPBC, are added for completeness. Setup names refer to the naming 
scheme used within the hotbit code. 24 ' 30 



Illustrations 



Examples for structures and distortions 



Bravais lattices (setup: Bravais) 

• conventional crystal lattices 

• translational symmetry in one or more dimensions 




Setup: Wedge 

• finite, symmetric molecules, clusters, marco- 
molecules, and viruses with cyclic axes (e.g. giant 
and high genus fullerenes, nanocones, dendrimers) 

• stretching and bending of one-dimensional struc- 
tures (nanowires, fibers, rings, tori, rod micelles, 
polymers, viruses) 

• bending of planar structures (ribbons, membranes, 
graphene, thin sheets) 

• adsorption and catalysis on cylindrical surface 




Setups: Chiral 



stretching, winding and unwinding of helical struc- 
tures (chiral and achiral nanotubes, springs, coils, 
helices, DNA, proteins, polymers, viruses) 
twisting of one-dimensional structures (tubes, wires, 
fibers, ribbons) 



Setup: DoubleChiral 

• calculating flat or twisted one-dimensional struc- 
tures with additional (generalized) reflection sym- 
metry 




Setups: Sphere and Saddle 

• approximate treatment only (small curvature) 

• spherical wrapping of planar structures (graphene, 
thin films, lipid membranes, vesicles and micelles) 

• distorting planar structures to a negative Gaussian 
curvature 

• adsorption and catalysis on distorted surface 



Setup: Slab 



one- and two-dimensional wires and slabs with (gen- 
eralized) reflection symmetry 
symmetric finite molecules and clusters 
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the smallest unit cell. Representative examples are car- 
bon nanotubes, where helical transformation have been 
discussed before 2 5 . The symmetry implies condition 
S^ 2 — S l i , and hence the restriction ft 2 M 2 = fti/i+27rra, 
where either n\ or n 2 is sampled freely. Important special 
case comes from \2 = tt, enabling, for example, an effi- 
cient simulation of twisted ribbons 17 . Yet another special 
case arises from T2 = 0, where helical transformation is 
extended by an independent rotation. 



E. Slab: surfaces and edges 

The transformations in this setup are two translations, 



S^r = r + mT 2 , (32) 

<S 2 n2 r = r + n 2 Ti, (33) 

and a reflection with partial translations, 

<S 3 n3 r = a n3 r + nn 3 Ti + r 2 n 3 T 2 , (34) 



where n, r 2 are half-integers and a n3 r ee [r — 2(r-z)z] ns is 
a reflection in x?/-plane. This setup is for 2D surfaces and 
ID edges (T 2 = 0), and allows halving the atom count. 
Values of (ti,t 2 ) are determined by lattice structure of 
the system in question, being either (0, 0), ( x /2, 0), (0, Y2), 

or (V2, 1 / 2 )- The y im Piy the identity <Sf = <S 2ri <S 2 T2 , 
leading to condition ft 3 = T\K\ + r 2 ft 2 + irm; k>i and 
n 2 can be sampled freely, while sampling of ^3 is given 
by m = 0, 1. Two values for ft 3 therefore introduce an 
additional component in two-dimensional band structure 
plots. 



F. Sphere: spherical symmetry 
(approximate treatment) 

As discussed in Ref. 20, symmetries can be treated 
also in an approximate fashion. Since for small angles 
rotations commute to the linear order, 7£(ai)7£(a 2 ) ~ 
7^(0:2)^(^1), two rotations around different axes but 
same origin, 

S^r = R(mai)r, (35) 
<S 2 n2 r = ft(n 2 a 2 )r, (36) 

approximately represents a spherical system with simu- 
lation cell of a square-conical shape. The ^-points can 
be freely sampled. This setup was recently used to cal- 
culate mean and Gaussian curvature moduli of graphene 
layers 20 . 



G. Saddle: negative Gaussian curvature 
(approximate treatment) 

Similarly, two rotations around different axes, 

<S x ni r = %ai)r, (37) 
<S 2 n2 r ee K(n 2 cx 2 ) (r - R ) + R , (38) 



such that the origin of first rotation (origin of the coor- 
dinate system) and the origin of second rotation, R , are 
located on different sides of the surface, approximately 
represent a saddle system. The ^-points can be freely 
sampled. This setup can be used to study membranes 
with negative Gaussian curvature in approximate fash- 
ion 20 . 



H. On choosing the symmetry setup 

Above all, various symmetry setups introduced by 
RPBC offer flexible access to different geometries, open- 
ing possibilities for distortion studies (see Table II). Be- 
sides, using symmetry setups for symmetric structures — 
whether classically or quantum-mechanically — bring re- 
duction of simulation costs; atom count reduction is pro- 
portional to the number of simulation cell images used in 
the setup. 

Because RPBC approach is centered around the sim- 
ple coordinate mapping r' = <S n r, computational over- 
head introduced by RPBC, as compared to translation- 
periodic implementations, should be in principle small; in 
a non-orthogonal tight-binding implementation hotbit, 
the computational overhead turned out negligible. 

It is important to choose the setup and the symmetry 
that fit the system best. On one hand, too symmet- 
ric setup will restrict natural conformations — regarding 
statistical sampling, for instance, RPBC faces the same 
artifacts as PBC, should the simulation cell be too small. 
On the other hand, exploiting too little symmetry will 
result in unnecessary simulation costs. 

Finally, general symmetries may require tighter geom- 
etry optimization criteria in case of small rotational an- 
gles. Consider, for example, Figure 2b and the total force 
on the atom in the simulation cell. Although the direct 
forces due to neighboring images may be large, the small 
wedge angle projects the combined force to have only a 
small radial component; geometrical arguments suggest 
that the radial component is diminished by a factor a 
(with a measured in radians), as compared to the direct 
forces. In other words, radial forces, although small, may 
cause certain type of center-of-mass creeping in radial di- 
rection, requiring more lengthy and demanding optimiza- 
tion process. With translation symmetry the center-of- 
mass movements never cost energy. 



V. TOWARDS PRACTICE: WARM-UP USING 
CLASSICAL POTENTIAL 

Let us first illustrate the practical aspects of the 
RPBC formalism with a classical treatment. Consider 
the interaction via pair potentials Vij(Rij), such that 
Vji(Rij) = Vij(Rij) an d Vjj ee for all / and J, where 
Rxj is the length of vector R/j = Rj — R/ separating 
particles / and J. The energy of the entire system, in- 
cluding all J\f images of the unit cell, with TV particles in 
each cell, is 

# = §XVjj(i*Jj), (39) 

IJ 
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where / and J still run over K = AfN particles. In order 
to reduce this expression to a single simulation cell, let 
/ and J run over N particles in one arbitrarily selected 
cell, while introducing an index n that runs over all its 
M images (J^ n 1 = A/"), giving 



E' 



N J\f 
I J nm 



Rn _ on/p m— n -ry \ — <?np m 
I — \ J — I J = O JrLj j 



(40) 



where 



Then 

Hjj = Rj — R/ is the vector separating the n th image of 
particle J and particle /. Since transformations <S n are 
isometric, the energy expression simplifies further to 



E' 



N JV 



I J nm 



N AT 



u 



Thus, energy per simulation cell is 



(41) 



I J n 



Here and below, indices / and J run over N particles in 
the selected simulation cell only; n runs over all unit cells 
where atom / at R/ still interacts with atom J at R™. 
Comparing energy expressions (40) and (41), one notices 
the obvious gain in the RPBC approach: Eq. (41) has 
one summation less. 

To calculate forces as derivatives of energy with respect 
to particles' positions, regard R n as a function of R/. In 
general, the partial derivatives of R n are elements of the 
rotation matrix 7£(u; n ) , 



[%)] 



(42) 



where i and j stand for the Cartesian components. Start- 
ing from Fj = — Vj-E, we first arrive at 



J, n 



J I 



I^Wlii^}, (43) 



where V is the derivative of V and R 7J is the unit vector 

Hj^/Rtj; note that R 7J ^ — Rjj. The second term in 
Eq. (43) may seem to be different from the first one, but 
it is not; this can be seen by observing the property 



R 



ji 



-^(u; n )R 7J 



(44) 



and by using the orthogonality of rotational matrix, 
[R-(wn)]ij = [TZ(u)—n)]ji to simplify the double sum in 

the second line of Eq. (43) to — R 7J . Thus, the total 
force acting on ion / in the simulation cell is 



(45) 



• 


• 


• 


• 


a 


- • 


• 




• 


b 




7f 





FIG. 2: (Color online) Peculiarities of symmetry setups il- 
lustrated on systems with one particle the in simulation cell 
(shaded), (a) While for translational symmetry the total force 
exerted on any atom is zero, (b) in general the total force (blue 
or light-gray) can be nonzero: a single particle in simulation 
cell can perform work on itself. 



This expression looks a simple as it should be: the total 
force exerted on atom / is the sum of the forces from all 
images of the other atoms J. 

Note that the expressions (41) and (45) are already suf- 
ficient to use RPBC in molecular dynamics with classical 
pair potentials. Generalization to three- or many-body 
terms is straightforward; a bond-order -type energy ex- 
pressions, for example, would yield in expressions like 



1 



N 



(46) 



UK nm 



VI. SOME PECULIARITIES OF REVISED PBC 

Differences between orientations of simulation cell im- 
ages lead to properties worth noting. Since we are 
all used to translational setups, these properties appear 
somewhat peculiar, although their origins are natural. 
Although applicable also in quantum mechanical meth- 
ods, for simplicity we illustrate these properties here us- 
ing classical pair potentials. 

First, consider the force exerted on the n th image of 
atom /. It can be shown to be symmetric in the sense 
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that F" = S n Fj [Note here that action of <S n on forces 
goes as for r i2 = r 2 — ri: the operation of <S n , loosely 
defined for both ri and r 2 by <S n r — 7£(u> n ) r + T n , gives 
<S n ri2 = 7^(u? n )ri2, as T n 's cancel out]. In the case of 
translational symmetry, the force exerted on given par- 
ticle in any cell is the same (lengths and directions of 
vectors F™ and Fj are the same). However, for symme- 
tries that contain rotational transformation, we have 

F? = ft(u> n )F 7 , (47) 

i.e. the vector of the force in n th cell is rotated by ou n 
with respect to the force in simulation cell. In wedge 
setup, for example, the forces resulting from cyclic po- 
tential demonstrate cyclic pattern. 

Related to this, consider a system with rotational ge- 
ometry and one atom in the simulation cell. Contrary 
to translations, when the forces exerted on any atom of 
a crystal with a monoatomic basis always cancel each 
other (sum up to zero; see Figure 2a), in the wedge setup 
it turns out that the resulting force exerted on the atom 
by its images may differ from zero; this is evident in Fig- 
ure 2b. Therefore, a single atom in the simulation cell 
can perform work on itself. In general, the force exerted 
on an atom due to its own images is 

n 

it may differ from zero if rotations or reflections are in- 
volved. 

Further, consider a system with wedge geometry and 
two atoms (/ and J) in the simulation cell. Now, the 
total force exerted on each of these atoms consists of two 
contributions, 

F/ = f/,j + f/,/, (49) 

and 

Fj = (j t i + fj,j, (50) 

where 

fv = Wj(* n ( 51 ) 

n 

is the force exerted on atom / due to atom J and all the 
images of J, and 

tj,i = J2 V rA R u)^n (52) 

n 

is the force exerted on atom J due to atom I and all the 
images of I [note that for a many-atom system, additional 
terms accounting for all the other atoms would appear in 
Eqs. (49) and (50)]. If only translations are involved, ro- 
tation matrix becomes the identity matrix, Eq. (44) gives 

R J7 = — R 7 j , and consequently f/ ? j = — fjj. It is 
tempting to regard this accidental property as Newton's 
third law. In general, however, we can have f j 5 j ^ — fj,j, 
but notice that this property is just an artifact of atom 
indexing. Namely, by atom / we mean atom / and all 
its periodic images, collectively. As defined by Eqs. (51)- 
(52), the forces f/ ? j and f j ; j are not the forces between 
two atoms / and J and, therefore, they do not have to 
follow Newton's third law in an atom-indexing sense. 



VII. REVISED PBC WITH A PRACTICAL 
QUANTUM METHOD 

Here we present the RPBC formulation with a 
practical quantum-mechanical method; we use a non- 
orthogonal tight-binding (NOTB) as our method of 
choice to illustrate the main concepts. To include 
the main aspects of the present-day standard density- 
functional (DFT) calculation, we discuss also charge self- 
consistency and use the concepts from self-consistent 
charge density-functional tight-binding (DFTB), which 
has been shown to be a rigorous approximation to full 
DFT. 25 As it will turn out, for an NOTB method RPBC 
merely introduces a simple transformation for the Hamil- 
tonian and overlap matrix elements. The electrostatic 
charge fluctuation correction is discussed separately in 
the next section. The NOTB-related classical repulsion 
energy has the form of a pair-potential, which was al- 
ready discussed above. But we start by giving short com- 
ments regarding practical basis sets. 



A. Comments on basis sets 

Numerically, there are many methods to solve the 
Schrodinger equation (2), and most of them can be used 
with RPBC — but some are more practical than others. 
Real-space grids and especially local basis are particu- 
larly suitable for RPBC implementation because they 
have more freedom to be rotated and adapted to vari- 
ous symmetries. This paper illustrates RPBC with local 
basis. On the other hand, plane waves are, by their very 
name, unsuitable for RPBC at least as such. 

The plane- wave method relies on building a basis from 
plane waves, which are constituted by the factor e* kr in 
Eq. (19). The wave function and potential V(r) are ex- 
panded in a Fourier series, in terms of ek(r) = e* kr , in or- 
der for the Schrodinger equation (2) to attain a practical 
algebraic form. 22 Yet Eq. (17) suggests that a similar type 
of method could be established also in the general case. 
Namely, one can build "symmetry-adapted plane waves" 
from the pre-factors e ^* n ( r ), with function n(r) chosen 
such that it properly reflects the symmetry of the system. 
Once the wave function and potential V(r) is expanded 
in a series with respect to these symmetry-adapted waves 
e w " n ^, however, the Schrodinger equation (2) not nec- 
essarily attain a neat algebraic form. Unraveling these 
issues are work-in-progress, and it remains a question 
whether this basis is numerically practical. On the other 
hand, for analytical investigations Eq. (17) is useful. 



B. Introducing notations: a short primer to NOTB 

The quantum-mechanical method that we use with 
RPBC is non-orthogonal tight-binding (NOTB) in gen- 
eral, and self-consistent charge density-functional tight- 
binding (DFTB) in particular. 26 30 For establishing no- 
tations, we hence review DFTB very briefly; readers fa- 
miliar with it can proceed directly to the next section. 

To represent a single-particle state ip a , the DFTB uses 
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a (typically minimal) local basis of ip^ (LCAO ansatz), 
V>a(r) = $>>M(r)> (53) 

where ^/i( r ) = ^^( r — IVr) is the orbital belonging to the 
atom I (fx e I); <^(r) is localized at the origin, oriented 

with respect to fixed Cartesian coordinates. 30 . The total 
energy is 

a (54) 

IJ KJ 

where f a is the occupation of the state ijj a and the Hamil- 
tonian matrix elements H^ u = ((p^\H[no] \ip u ). (Notation 

H, instead of the usual H°, is chosen to avoid superscript 
confusions later. 25 ) In the second term the excess Mul- 
liken charges Aqi = qj — q®, as compared to charges q® 
of valence electrons from neutral system, interact with 
effective electrostatic interaction 7/j. 

The efficiency of the DFTB relies on the fact that ma- 
trix elements H^, 

= Jd'r^ir-R^H^ir-Rj) (55) 

with \± G /, v G J, as well as S , are tabulated with re- 
spect to distances between orbital centers Rjj = |Rj — 
R/|. The matrix elements with all the possible orienta- 
tions of the orbitals <p°* and (p® are accounted for by a 
superposition of the matrix elements with Slater-Koster 
transformation tables 31 and just a few basic matrix ele- 
ments. 



Single-particle states 

^(r) = ]Tc»^ K (r) (59) 

/I 

then fulfill the revised B loch's theorem, Eq. (6), by con- 
struction. 

The image orbital <^"(r) is the simulation cell orbital 
<£/x(r) transformed by £)(<S n ), and is not an eigenfunction 

of D(S n ). When acting on the orbital ^(r), located on 
atom / at R/ and oriented along Cartesian axes, transfor- 
mation D(S n ) can be thought to be the superposition of 
three consecutive transformations of i) translation of the 
orbital from atom / to origin [becoming hence the orbital 
(/^(r)], ii) rotation of the orbital though a n -radians, to 
adapt to the "orientation" of the image n, and iii) trans- 
lation of the orbital to R" , the position of atom / in the 
image n (Figure 3a), i.e., 

D(S n ) = D(T n ")D(n(cx n ) )I)\T^ A ). (60) 
This gives 

^(r) = J2 ^(ft(«n) ) ¥V ( r " ( R " " *V))> (61) 

where yy(r - (R°, - R M >)) = <^°,(r - R"), with R^ 
meaning R/ such that \x G /, and 

£>°„(^(a n ) ) = <$|£(ft(a n ) )|^> (62) 

being the representation of rotation D(TZ(ot n ) ) in basis of 
orbitals ^°x( r ) located at the origin with the conventional 
orientation along Cartesian axes; the explicit form for 
D^ u (7l(cx n ) ) is given in Appendix A. The orientation of 
the orbitals ^™( r ) i n different images n hence follow the 
symmetry, as illustrated in Figure 3b. 



C. Revised Bloch waves 

The Bloch basis functions for translational symmetry, 
or Bloch waves, are given with a local basis ^(r) by the 
familiar expression 30 

T n 

where T n (p^(r) = (/^(r — T n ); Bloch waves extend 
through the entire system. Analogously, for general sym- 
metries, atomic orbitals of all M images, 

^(r) == D(S> M (r), (57) 

are assembled to create new basis functions according to 

^«( r ) = VN E eiK ' n D(S^,(v), (58) 

n 

with <^«(r) = (r|^). 



D. Transformation of matrix elements 

In the revised Bloch basis, Eq. (58), the Hamiltonian 
and overlap matrix elements H^^k^k') = ((p"\H\(p* ) 
and S^ v {n^ k') = (<^|<^£ ) transform as well; we need to 
account for orbital rotations. 

While here we consider Hamiltonian matrix elements 
only, similar expressions hold for overlap matrix elements 
as well. Using the explicit form (58) for new basis func- 
tions, we obtain 

H» v {K,K') = W\H\rt') 

= ^VJ e ---+-'-"^(m,n), 

m,n 

where 

^(m,n) = <^|tf|^> (63) 
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FIG. 3: (Color online) (a) Illustration of Eq. (60): three- 
step transformation of orbital ^(r) to orbital <^"(r). (b) 
Illustration of orbitals <p"(r) for the entire system in wedge 
setup. Orbitals are symmetrically aligned (simulation cell is 
shaded) . 



(H and H transform the same way, so we use H instead 
of H). Since H and D(S n ) commute, (<p™\H\<p?) = 
(ip \H\p™~ m ), or with a shorter notation 

i^„(m, n) #^(0, n - m) = H^ v {n - m). (64) 

Hence the matrix elements become diagonal in as ex- 
pected, 



^( K ) = ^e iK - n ^(n), 



(65) 



and the ^-dependence of matrix elements ends up resem- 
bling translation-periodic formalism 30 precisely — only k 
has changed to k. 

However, the matrix element, Eq. (64), is now given 

by 



n) = |d 3 r^( 



(r)H D(S a )Mr) 



(66) 



and becomes hence more complex due to rotations in- 
duced by D(S n ). By substituting Eq. (61) inside the 
square brackets, a short manipulation yields 



(67) 



where 



H: 



;» = J d 3 r ^(r)iJ^(r - (R™ - R„)) 

= yd 3 r^(r-R M )H^(r-R^) (68) 

is the matrix element involving translations only, with 
orbitals (p^ , centered at R M , orbitals ip° , centered at the 
origin, all orbitals having a fixed orientation. (Orbitals 
oriented with respect to the fixed Cartesian coordinates 
of the simulation cell.) The matrix element H'^ v as given 
by Eq. (68) is practically the same as the one in Eq. (55), 
and it is straightforward to calculate in any tight-binding 
implementation. Eq. (67) is a simple matrix multiplica- 
tion, and in compact matrix notation 



H(n)=H'(n).D°(ft(a n )). 



(69) 



Thus, compared to the translation-periodic formalism 30 , 
matrix elements need simply to be transformed by 
D°(7£(a n ) ). Eq. (69) is valid within the two-center ap- 
proximation, otherwise one needs to return to Eq. (66). 



E. Total energy and forces in NOTB 

The total energy expression for NOTB with revised 
PBC is therefore 



£tot=J>" Tr [/>(*)■ H(tt) 



KJ 



(70) 



where w n are weights for equivalent ^-points (^2 K w n = 
1), and 

P ^(n) = J2fa(n)[c a ;(nK(n)] (71) 



is the ^-dependent density matrix. Mulliken charge anal- 
ysis is done as before (although with transformed S- 
matrix), and the electrostatic charge- fluctuation correc- 
tion remains the same, 



E 2 nd = lJ2 GljAqiAqj 



IJ 



except that the expression 

n 

is different, now using the generalized symmetry. 



(72) 



(73) 
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Minimizing E tot -^2 a f a (K)e a (n) [(^ a \^a) ~ 1] with re- 
spect to c a * (k) yields the generalized eigenvalue problem 

^c^(K)[H^(K)-e a (n)S lxu (n)] = 0, VK,/x,a, (74) 



where 

)J2U G lK + G JK )Aq K 
K 

for ji G I and v G J. The secular Eq. (74) looks very 
familiar — RPBC does not alter the appearance of the 
NOTB formalism at all. 

Differentiation of the band structure energy with re- 
spect to the atoms' positions yields band-structure forces 
straightforwardly as 

F/ = - w * Tr ' [ dH («M«) - dS(K)p en («)] . (75) 



Here Trj stands for partial trace over orbitals of atom / 
only, 

/C» =^2fa(K)e a (K) [<£*(*)<£(«)] (76) 



is the energy-weighted density matrix, and 

dH^(K) = - e iKn Vj^(n), ueJ (77) 

n 

with the same for dS(tt). Gradients Vji^ /izy (n), 
Vj5 /Ul/ (n) can be straightforwardly calculated from the 
gradients of H'^(n) and 5^(n), Eqs. (67)-(68). 

We still wish to emphasize the obvious result: while 
internal structure of the Hamiltonian and overlap matrix 
elements H^ u (n) and S^ u (k) for general symmetries are 
different from the translational case, the final equations 
are trivially the same — k only changes to k. 

VIII. ELECTROSTATICS 

A crucial component of the interatomic interaction in 
ionic or partially polarized solids is electrostatics. Naive 
direct summation of the long-ranged electrostatic pair 
interaction, while possible using the approach outlined 
in Sec. V, is slow and only conditionally convergent and 
hence numerically prohibited. 

Methods with better convergence properties have been 
established for periodic systems with translational sym- 
metry in three dimensions. Most of these are based on 
Ewald summation 32 which do not straightforwardly gen- 
eralize to partially periodic systems. Hence, we use a 
direct summation method whose convergence is accel- 
erated by the use of a hierarchically telescoped mul- 
tipole expansion 33 , in the spirit of the Fast Multipole 
Method (FMM) 34 36 . In a multipole expansion scheme, 
the charge distribution p and electrostatic potential $ 
are expanded around suitably chosen origins. These ex- 
pansions are called the multipole expansion M and the 



multipole-to-local 



a = i i | I I M I i i n 

— mu|t 'P o|e " to_mu|t 'P o|e 

a = 1 j I I j I I j I I j I I I I | I I | I I | i n 

a = 2 | | | 1 

a = 3 j j 



FIG. 4: (Color online) Sketch of the operations needed for 
the telescoped fast-multipole expansion demonstrated for an 
example one-dimensional Bravais lattice and Ni = 3. The 
value of a denotes the level of the expansion. The method 
computes the electrostatic potential and field within the sim- 
ulation cell (yellow), multipole-to- multipole operation: 
At a — 0, the multipole moment Mo is the moment of the 
simulation cell (yellow). At each successive level M Q is the 
combined multipole moment of N (here N = 3, red) units 
of level a — 1, where the basic repeat unit at level a is de- 
noted by the broken lines, multipole-to-local operation: 
At each level a, the local expansion at the origin is deter- 
mined from the multipole expansion at level a (from the green 
cells), local-to- local operation: Within the simulation cell 
(yellow), the local expansion is shifted to each particle to de- 
termine the electrostatic potential and field at the particle's 
position. This is then combined with the near-field contribu- 
tion (from the gray repeat units). 



local expansion L, respectively. We will carry out the ex- 
pansion in terms of spherical harmonics, but also Taylor 
expansions have been reported 37 . 

The determination of the electrostatic potential and 
field requires three operations (see also Figure 4). First 
operation is to shift the origin of the multipole expan- 
sion; it is commonly called the multipole-to-multipole 
transformation. This operation is required to telescope 
the multipoles, i.e. to join a number of adjacent multi- 
pole expansions into a single one. Second operation is to 
determine the local expansion from the multipole expan- 
sion, the multipole-to-local transformation. This opera- 
tion solves the electrostatic problem and is typically the 
most computationally demanding step. The third and fi- 
nal operation is to shift the origin of the local expansion 
to arbitrary points in space, the local-to-local transfor- 
mation. 

Since the method operates in real space, RPBCs are 
straightforwardly implemented in such a scheme, because 
it requires only the mapping r' = <S n r. Basically, in addi- 
tion to translation (i.e. the multipole-to-multipole trans- 
formation) the multipole moments of the expansion need 
to be rotated. For spherical harmonics, such a rotation 
is described by the Wigner d-matrix 38 and hence the ro- 
tation operation is similar to the rotation of the orbitals 
described in section VII. The general scheme for electro- 
statics in RPBCs, that includes reflection and inflection 
operations, is described in more detail in the following. 
Figure 4 shows a graphical illustration of the individual 
operations. 

We would like to note that such multipole expansion 
scheme is also straightforwardly transferable to full DFT 
with Gaussian or other atom-centered basis functions. 
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An example of such a method for nonperiodic systems is 
the Gaussian very fast multipole method (GvFMM) 39 ' 40 
and the continuous fast multipole method (CFMM) 41 . 
Since in the far-field the only relevant moment of the 
Gaussian charge densities is the monopole, the far-field 
solution for Gaussians is identical to that obtained in the 
FMM for point charges. The near field solution contains 
two-electron repulsion integrals which need to be evalu- 
ated directly 42 . Indeed, the 7/j Coulomb integral of the 
self-consistent charge NOTB formalism is the two-center 
expression of such an overlap contribution for the inter- 
action of Gaussian charge distributions with s-symmetry. 



A. Multipole expansion 



The multipole moments are given by 



(78) 



where / denotes atoms and AR/ is the position vector of 
atom / computed with respect to some expansion origin 
R , i.e. AR/ = R/-R . Note that \m\ < I. The i?p[r) 
in Eq. (78) are the regular solid harmonics, 



i?r(r) = 



1 



(Z + m)\ 



r l e - im *P™ (cos 6), 



(79) 



Multipole-to-multipole, multipole-to-local, and 
local-to-local transformations 



The general transformations underlying the fast- 
multipole expansion method have been described in de- 
tail elsewhere 35 ' 47 and can indeed be found in textbooks 
on electricity and magnetism 43 . We will therefore only 
summarize the basic equations, and indicate where the 
symmetry transformation operators <S n need to be ap- 
plied. 

Multipole-to-multipole. Let Ro denote the origin of the 
multipole expansion M. The expansion M is located at 
0. Then, the expansion coefficients M™ are given by the 
convolution 



1 



M, m = EE M£iC7(R ) = M o R(R ) 

A=0 //=- 



(82) 



where M/ 71 = if \m\ > I. Eq. (82) defines the convolu- 
tion operator o. 

Direct evaluation of the telescoped sum of repeating 
cells requires to compute the joint multipole moments. 
Let M a be the multipole moment at stage a of that 
telescoping process (with M being the bare multipole 
moment of the simulation cell). The multipole moment 
at stage a is then given by the recursion formula 



^M a _ lOj D(5 N - n )R(0), 



(83) 



where P{ n {x) are the associated Legendre polynomials 
and r, 6 and cp are the length, inclination and azimuth 
angle of r, respectively. The corresponding conjugates 
are the irregular solid harmonics ij m (r) given by 

/ f m (r) = (Z - m)!r- ( ' +1) e im ^ m (cos 0). (80) 

The expansion of the electrostatic interaction into mul- 
tipoles is based on the identity 43 



00 A 



= E E ^rw W) 



r — 1" 



(81) 



A=0 jjb=-\ 



for the free-space Green's functions. The overall scheme 
is hence a real-space method. A multiplicative decompo- 
sition of the kernel such as Eq. (81) can also be achieved 
using wavelets 44 ' 45 that are indeed related to the fast 
multipole method 46 . It is hence conceivable that wavelet- 
based algorithms could be used here instead of the cur- 
rently proposed scheme. 

Note that the normalization used in Eqs. (79) and (80) 
is non-standard. The usual solid harmonics are given by 
yjjl — mjUj + m)\Ry i (r). This particular choice of nor- 
malization is motivated by numerical considerations. At 
large expansion orders I the value of r l becomes large and 
the factor (Z + m)\ provides some compensation. While 
this is not a rigorous approach towards numerical stabil- 
ity, in all practical situations the expansion is cut-off at 
a value Z max . The approach outlined here is stable up to 
at least Z max = 10. Well-converged solutions are typically 
obtained for Z max = 8. 



where n = (711,712, . . .) is the usual symmetry index and 
the sum runs from —(N{ — l)/2 to (N{ — l)/2 for each 
rii. Furthermore, N a n denotes a component- wise multi- 
plication and N a = ((Ni) a , (N2) 01 , . . .), hence summing 
up the contributions of N{ repeat units for the symme- 
try Si. Note that N a gives the distance of neighboring 
repeat units at stage a (in numbers of simulation cells). 

Since D is an affine transformation we can split it ac- 
cording to D = T + Dq [see also Eq. (25)], where T con- 
tains the translation operation and Dq the linear trans- 
formation. We now absorb the translational contribution 
into the regular solid harmonics and rotate the multipole 
expansion. Eq. (83) then becomes 



Nq 



')M a _i oR(5 ° n O) 



(84) 



where <S -Nc * n is the Cartesian distance between indi- 
vidual repeat units at stage a. Intuitively, the multipole 
moment of a single stage in the telescoped expansion has 
to be given by a combination of rotated moments, which 
leads directly to Eq. (84). 

Multipole-to-local. Let the multipole expansion M be 
localized at 0. Then the expansion L of the local poten- 
tial at Ro is given by 

LT = Y.T. M A W(Ro) = [M * Wr , (85) 

A=0 /x=-A 

where * is another convolution operator. In analogy to 
Eq. (84) the local expansion of the potential at the center 
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of the unit cell is given by 

« max -2 

L r= E X)^o( 5 " Non ) M «* I (' s " N,,no )' ( 86 ) 

a>=0 n 

where the sum over n now runs from —Ni to N{ excluding 
the terms where the well-separateness criterion does not 
hold, i.e. where |n^| < Ni ~ 1 for any i. In total, this sums 
up 7V" max images of the simulation cell for symmetry Si. 

Local-to-local. Let the local expansion L be located at 
0. Then the expansion L located at Ro is given by 

LT = E E L^YR^Ro). (87) 

A=0 fi= — X 

Note that since the local expansion is only applied within 
the simulation cell, no symmetry operation needs to be 
applied here. 

Near-field. The contribution of the neighboring Ni re- 
peat units to the electrostatic potential is computed by 
direct summation, as discussed in detail in Sec. V. 



C. Transformation of multipole moments 

We apply a general linear transformation C to all coor- 
dinates. In particular, in 7Z 3 such linear transformation 
has a representation in a 3 x 3 matrix T and hence the 
transformed positions rf are given by rf = £r^ = Tr^. 
Under this transformation, the solid harmonics transform 
as ^{OjWp^i) — ^{TYi) and from group-theoretical 
arguments 21 it becomes clear that C has a block diagonal 
matrix representation in the function space of the solid 
harmonics. In particular, 

i 

D(C)Rr(r i )= J2 D^ m ,Rf(r i ) = m(T)R l (r i ) 

m' = -l 

(88) 

for each I and m with \m\ < I. (Compare also Eq. (69) 
for the rotation of Hamiltonian and overlap matrix ele- 
ments.) The computation of matrix up to arbitrary 
order in I is described in Appendix B. 



IX. APPLICATION EXAMPLE: POLYALANINE 

We are now in a position to use the RPBC in combina- 
tion with the self-consistent charge tight-binding model, 
using polyalanine in vacuum as an example system. This 
chiral molecule belongs to the class of polypeptides, 
where the peptide unit alanine has a single methyl side- 
group. The RPBCs allow to continuously search for con- 
formations of this molecule without the need to choose 
conformations that are commensurate with a Bravais unit 
cell 48 . 

The molecule is modeled using chiral symmetry. We 
use the O-C-N-H tight-binding parametrization of Elst- 
ner et a/. 49 , augmented by van-der-Waals interactions 50 
using the polarization data of Miller 51 , and sample n- 
space using 20 equally spaced ^-points. The electrostat- 
ics interaction is computed using / max = 8, N\ = 11 and 




FIG. 5: (Color online) Conformations of polyalanine in vac- 
uum, (a) 7r-helix and (b) a-helix. Atoms colored in black 
show the simulation cell consisting of 10 atoms. Blue dot- 
ted lines indicate hydrogen bonds between the oxygen and 
hydrogen atoms. 



c^max = 5. For each conformation, we vary the recurrence 
length along the molecular axis z and optimize the chiral 
angle x- 

The equilibrium conformations obtained for the 7r-helix 
and the a- helix are shown in Fig. 5, where the atoms 
shown in black highlight the atoms within a single unit 
cell. We find an equilibrium repeat length along the 
molecule's axis of 1.21 A and 1.57 A for the n- and a- 
helix, respectively, which corresponds to a chiral angle of 
82.5° and 100.84°. This compares well to full density- 
functional calculations where a repeat length of 1.17 A 
and 1.50 A is found for the n- and a- helix, respectively 48 . 

Note that using a Bravais unit cell the chiral angle is 
given by = 360°m/7V, where m is the number of chain 
turns and N the number of residues per unit cell. Hence, 
for m = 1 the peptide's chiral angle is fixed to either 
90° or 120°. If the goal is to bracket the chiral angle of 
the a-helix in an interval smaller than 1° the Bravais cell 
needs to contain at least m = 32 turns and N = 114 
residues, compared to a single residue when using the 
RPBC approach. 



X. CONCLUSIONS 

We have demonstrated that RPBC is a simple but pow- 
erful atomistic simulation technique with easy implemen- 
tation. It brings considerable efficiency gains and creates 
an access to distorted material simulations. Implementa- 
tions of RPBC in quantum and classical codes will open 
possibilities for new types of simulations, such as stud- 
ies on crack propagation under complex load conditions, 
DNA stretching, catalysis on curved surfaces, and stud- 
ies of macromolecules in nanochemical and nanobiolog- 
ical systems. We acknowledge that an implementation 
of RPBC in plane- wave DFT codes is not straightfor- 
ward and requires further theoretical research. Never- 
theless, the present paper shows that an implementation 
in self-consistent charge NOTB is possible and it contains 
all the ingredients for an implementation into full DFT 
codes that employ local basis sets and the fast multipole 
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method, such as for example Q-Chem. 52 We believe, that 
RPBCs, both in classical and quantum- mechanical for- 
mulation, has the potential to become a truly useful tool 
in computational nanoscience. 
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Jyvaskyla for hospitality during a visit that was funded 
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Appendix A: The transformation matrix D°(7^(a n )) 



where Dq are (2Z + 1) x (2Z + 1) matrix blocks that rotate 
orbitals with angular momenta I = 0,1, 2,.... Since the 
orientation of s-orbitals does not change upon rotation, 

= 1. Since tight-binding basis orbitals are chosen 
real (see Ref. 30), three p-orbitals, p x -,PyiPzi transform 

as x, y, z coordinates: Dq 1 ^ = 1Z(cx n ) . 

One way to obtain general expression for D°(7£(a n ) ) 
in tight-binding basis (Al) is to carry out similarity 
transformation of the rotation representation, which in 
commonly used basis of angular momentum eigenstates 
|Z,ra) is given by 23 ' 38 



a£U«,&7) = e- im ' a d%, m {P)e- im \ (A3) 



In quantum mechanics, transformation of wave func- 
tions {e.g. orbitals) upon rotation is found in group rep- 
resentation theory. For atom's basis orbitals ordered as 

{K)} = {\S), \Px),\Py),\Pz), 

), \d xy ), \d yz ), \d zx ), \d z 2), . . . }, (Al) 
rotation representation has block-diagonal form 

/ 



D°(ft(a„)) 











...\ 





D (i) 













D (2) 










' • / 



(A2) 



where d 



(0 

rn'rn 



are Wigner d-matrices and a, /?, 7 are Eu- 

(2) 

ler's angles. For d-orbitals, explicit expression for Dq ; 
becomes tedious; we omit it here. 

For illustration, however, we present a special case rel- 
evant for wedge and chiral setups. For rotations about 
z-axis, p-orbitals transform with 



cos(a n ) — sin(a n ) \ 
sin(a n ) cos(a n ) , (A4) 
1/ 



and d-orbitals transform with 



cos(2a n ) - sin(2a n ) \ 

sin(2a n ) cos(2a n ) 0* 

cos(a n ) sin(a n ) 

— sin(a n ) cos(a n ) , 

1/ 



(A5) 



note that it is the transpose of these matrices that is 
acting in Eq. (61). 



Appendix B: The transformation matrix D^(T) 



note the difference between D^' and D^, Dq from pre- 
vious section. 

The matrix for I = and I = 1 can be evaluated 
explicitly. In particular Z)q q = 1 and 



The matrices D^) are related to the Wigner d-matrix, 
which describes the rotation of spherical harmonics, with 
normalization factors: 



D 



(0 _ Kl-m')\(l + m')\ (Q . 

" V (l-m)l(l + m)l mm " [ } 
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1 / Tyy ~\~ T XX -\~ i(T X y Tyx} T XZ iTy Z Tyy T XX ~\~ %(T X y ~\~ Ty #)\ 

D (1) = - ( 2(T zx + iT zy ) 2T ZZ 2(-T zx + iT zy ) 1, (B2) 

\ T XX ~\- Tyy i(T X y ~\~ Ty ^ ) T X% ITy Z T XX H~ Tyy l(T X y Tyn^) J 

where the latter one is the equivalent of Eq. (A4) for the transformation of p-orbitals around z-axis. 

In order to evaluate the matrices for I > 2 up to arbitrary order in Z, we use a recursion scheme that has been 
developed for the rotation of spherical harmonics 53 ' 54 . In terms of the regular solid harmonics given by Eq. (79) the 
recursion becomes 

f)(i) ( -i.n+ifS(i) n^- 1 ) +o ' n D (1) n^- 1 ) +o 1 ' n - 1 D (1) b {l - 1] \ (BD 

a l-i V 7 

_/ <m< /.p -^f.-WnW n^- 1 ) , 0,n n(l)n(I-l) l,n-l^(l)^(Z-l) \ 

I < 171 < I . U mn — Q I a i _ 1 ^0,-1^771,71+1 + a l-l- U 0,0- U m,n + a £-l ^0, l^m, Ti-1 J l 154 ^ 

f)(0 _ 1 n(!) n^- 1 ) i w o,n f>(i) f>(z-i) , i,n-i ~(i) ~(z-i) \ 

- i /_! l a Z-l ^1-1^-1^+1+^-1^1,0^-1,71+^-1 U \,\ U l-\,n-\) -> 
a l-l V 7 



where, as opposed to a rotation of spherical harmon- Equations (B2) to (B8) are the explicit representation of 
ics 53 ' 54 , Eqs. (B3)-(B5) hold for general linear transfer- the rotation operation. 



mations C with det T ^ 1 (i.e. inversion or shear). The 

a 
I 



prefactors af ,rn are given by 



a -i, ra Jl-m + W-m+l) (B6) 
o,m v T" t j-a 1 ~~ T" j-y 

a i = 07 I 1 ' ' " > 

- 2(2/+ i) • < E8 » 





2(2/ + 1) 








HI) 




21 + 1 






hm + 2)(7 + mH 


HI) 
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